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Background: Mean-field methods based on an energy density functional (EDF) are powerful tools used to 
describe many properties of nuclei in the entirety of the nuclear chart. The accuracy required on energies for 
nuclear physics and astrophysics applications is of the order of 500 keV and much effort is undertaken to build 
EDFs that meet this requirement. 

Purpose: The mean-field calculations have to be accurate enough in order to preserve the accuracy of the EDF. 
We study this numerical accuracy in detail for a specific numerical choice of representation for the mean-field 
equations that can accommodate any kind of symmetry breaking. 

Method: The method that we use is a particular implementation of 3-diniensional mesh calculations. Its nu¬ 
merical accuracy is governed by three main factors: the size of the box in which the nucleus is confined, the way 
numerical derivatives are calculated and the distance between the points on the mesh. 

Results: We have examined the dependence of the results on these three factors for spherical doubly-magic 
nuclei, neutron-rich ®^Ne, the fission barrier of ^"^^Pu and isotopic chains around Z = 50. 

Conclusions: Mesh calculations offer the user extensive control over the numerical accuracy of the solution 
scheme. By making appropriate choices for the numerical scheme the achievable accuracy is well below the model 
uncertainties of mean-field methods. 


I. INTRODUCTION 

The self-consistent mean-field approach, based on an 
energy density functional (EDF), is a tool of choice to 
study nuclei in any region of the nuclear chart [ij. It al¬ 
lows to calculate the properties of the ground state but 
also of alternative configurations, like shape isomers, or 
to follow the behaviur of a nucleus along rotational bands 
or along fission paths. Often, one is not directly inter¬ 
ested in the total binding energy of a specific nucleus 
but by its evolution along a series of isotopes or isotones, 
which can signal structural changes for given neutron or 
proton numbers. 

Motivated by the needs of the nuclear physics and as¬ 
trophysics communities, large efforts are underway to 
push the predictive power of nuclear mass models well 
below the 500 keV level. To reach this goal, the proto¬ 
cols used to adjust the EDF’s parameters have been re¬ 
visited. In particular, methods are being developed M 
to quantify the statistical uncertainty of these parame¬ 
ters. However, besides the errors on observables due to 
these uncertainties, there is also a numerical error due to 
the way the SOME equations are solved. One needs to 
verify that the numerics does not introduce errors that 
are larger than the maximum error tolerated for mass 
models. More importantly, these errors should not vary 
too rapidly from one nucleus to the other, to avoid a 
spurious behavior of mass differences. 
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The numerical methods used to solve the mean-field 
equations can be classified according to the way the 
single-particle wave functions are represented: by coor¬ 
dinate space techniques or by a basis expansion. Co¬ 
ordinate space techniques represent the single-particle 
wave functions in a discretized, finite volume. Several 
discretization techniques exist, utilizing finite difference 
formulas 13, Fourier transformations [Q, B-spline s H , 
wavelets 0 and the Lagrange Mesh method [lOl - H^ . 

The second family of numerical representations in¬ 
volves expanding the single-particle wave functions on 
some chosen (finite) set of basis states. Usually these 
basis states are harmonic oscillator (HO) eigenstates, al¬ 
though the details often vary. 

While the origin of numerical errors is quite different 
for both families of representations, the type of EDF does 
not seem to influence the accuracy of the methods much. 
The three main families (relativistic EDFs, zero-range 
Skyrme EDFs or finite-range EDFs) require similar num¬ 
bers of basis states to achieve a similar precision (see 
e.g. [3-[l3])- In what follows, we will limit ourselves to 
the study of zero-range Skyrme EDFs. 

It is the aim of this paper to study the numerical ac¬ 
curacy of a specific implementation of coordinate space 
techniques: representation on a three-dimensional carte¬ 
sian mesh of equidistant points. We will focus on two 
specific techniques: finite-difference (ED) formulas and 
the Lagrange Mesh (LM) method, which are the ones im¬ 
plemented in our codes. As far as we can infer from the 
tests published in the literature, the accuracy obtained 
with the other techniques mentioned above is similar to 
the one obtained within our LM scheme. Most of the in- 
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formation relative to the tools that we have developed has 
been presented for the particular implementation made 
in the code Ev8 More involved implementations 

have also been used, which differ from Ev8 only in that 
they impose less symmetries on the nucleus. The pres¬ 
ence of these symmetries in general allows for the reduc¬ 
tion of the dimension of the problem, e.g. in EV8 it allows 
for the reduction of the space by a factor of 1/8. 

The article is organized as follows: first we define pre¬ 
cisely the quantities that will be used to characterize the 
accuracy of a mean-field calculation. Next, we review 
the basic ingredients needed to define wave functions on 
a cartesian mesh and to calculate derivatives and inte¬ 
grals in this representation. We then discuss the main 
sources of numerical errors: the size of the box in which 
the nucleus is confined and the step size of the mesh. We 
discuss the numerical accuracy that can be achieved by 
comparing energies and radii of doubly magic nuclei with 
those obtained with a spherical code. Finally, we check 
the convergence of energies, radii and the multipole mo¬ 
ments of deformed nuclei by comparing results obtained 
with decreasing mesh discretization lengths. 


II. DEFINITION OF USEFUL QUANTITIES 

A mean-field configuration is characterized by its en¬ 
ergy, its rms radius and by multipole moments. In this 
section we define these quantities whose dependence on 
the mesh parameters will be studied. 


with coupling constants as defined in Ref. [l^. The ki¬ 
netic energy just depends on the kinetic density 


I?kin 



( 3 ) 


of protons and neutrons. While the Skyrme and kinetic 
energies are local functionals of the densities, the direct 
Coulomb energy is a nonlocal functional of the proton 
density Pp{r) 




Pp{r)pp{r') 
|r - r'l 


( 4 ) 


Compared to the other terms contributing to the total en¬ 
ergy o. the exact calculation of the Coulomb exchange 
energy is orders of magnitude more costly as it is a func¬ 
tional of the complete nonlocal one-body density matrix. 
As a consequence, the local Slater approximation that is 
of similar numerical cost as the Skyrme energy is used 
instead 

= (5) 

The pairing energy contribution to the energy is: 

^'pair = 'y ] fk UkVk fm UmVm 
k,m>0 


where the are anti-symmetrized matrix elements 

of the pairing interaction and the fi are cutoff factors, 
both of which are specified in Appendix [B1 

The expression for the c.m. correction, which is not 
relevant for our discussion, can be found in Ref. M- 


A. Total Energy 


B. Dimensionless multipole moments 


For a time-reversal invariant system as assumed here, 
the total energy is composed of the kinetic energy, the 
Skyrme energy describing the strong interaction in the 
particle-hole channel, the pairing ener gy, the Coulomb 
energy and a center-of-mass correction [l9| 

^tot — A^kin T Askyi-jne T Apair T AqquI T Acm ■ (1) 

For the parameterizations used throughout this article, 
the Skyrme EDF takes the form of the sum of various bi¬ 
linear combinations of the isoscalar {t = 0) and isovector 
(t = 1) local densities Pti^^), kinetic densities Tt(r) and 
spin-current densities Jt,pv(j)^ p, v = x, y, z, 


Askyrme — Ap2 Ep2+a Ep-j- -j- Ep/^p -\- Ep^j -j- Ejj 
= (cfbo] Pt + cfpo Pt + Cl Pt Tt 

4- _r» 1 


i=0,l 

"fZ 


+C^^Pt Apt + C'^-'^ptV-Ji 




t 

j 5 


( 2 ) 




As in [l9j |. the dimensionless multipole moments 
are related to the matrix elements of the multipole oper¬ 
ators Q(^ = Yfm(r) by 




Ltt 

SAp" 


{Qim) j 


( 7 ) 


where Aq = 1.2 fm. When m is omitted we imply it 
to be zero. 


C. Radii 


Another set of observables, related to the density pro¬ 
file of the nucleus, are the mean-square (ms) radii, root- 
mean-square (rms) radii and the isotopic shifts. The ms 
radius of the proton (g = p), neutron (g = n), and total 
density distribution is defined as 



( 8 ) 


d^r [p„(r) -Epp(r)] . 


(9) 
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The root-mean-square (rms) radii are then the square 
root of the corresponding mean-square radius. 

Similarly, we will present results for the isotope shift of 
charge radii that are calculated as the difference between 
the proton ms radius of an isotope with N neutrons and 
a reference isotope with Nq neutrons 

Sr^{N,Z)=rliN,Z)-rliNo.Z) ( 10 ) 

without any corrections. 


III. COORDINATE SPACE REPRESENTATION 

Assuming a 3-dimensional cartesian mesh, a function 
<l)(r) = $(a::, y, z) is represented by the tensor of its 
values at the collocation points {xp,yq,Zs) 

$(r) = {$(xp, 2 /g,z^)} = (11) 

A mesh can be defined in several ways, depending on the 
choice of the collocation points. For example, the origin 
of the coordinate system and the boundaries of the box 
can be included as collocation points or not. Different 
choices can also be made for the boundary conditions at 
the edge of the box. 

To set up the self-consistent mean-field equations, one 
has to vary the EDF with respect to This requires 
to define prescriptions to calculate derivatives and inte¬ 
grals from the values of on the mesh. Several choices 
for derivatives have been explored over the years. 


A. Derivatives on a mesh 

The most straightforward possibility to set-up a 
coordinate-space representation of the self-consistent 
mean-field equations is provided by the finite-difference 
method, a widely-used tool to solve partial differential 
equations [^ . 

In such a scheme, the derivatives are calculated with n- 
point finite-difference formulas, and the integrals are ob¬ 
tained by summing up the integrand at the mesh points 
multiplied by a suitable volume element. 

There are three factors that determine the accuracy 
that can be achieved with the finite-difference method. 
One is the overall resolution scale provided by the mesh 
spacing; decreasing the distance between mesh points im¬ 
proves the accuracy. Second, the higher the order of the 
finite-difference formulas used for a given mesh spacing 
the better the accuracy. In both cases, however, bet¬ 
ter accuracy means also increase of the numerical cost. 
Third, there are internal inconsistencies introduced by 
the method itself. For example, taking twice the nu¬ 
merical first-order derivatives of a given function is not 
equivalent to applying the numerical second-order deriva¬ 
tives. Also, the numerical derivatives are not the inverse 
of the numerical integration. It is only for very small step 
sizes well below 0.1 fm that these internal inconsistencies 


become irrelevant. While such small step sizes can be 
easily handled in spherical Id codes [2l|, the required 
storage is prohibitive in axial 2d and cartesian 3d codes. 
In addition, such step sizes are much smaller than what 
can be expected to be the physically relevant resolution 
scale, see for example the arguments brought forward in 
Ref. In. 

Several other schemes have been developed in the past 
with a better consistency between derivation and integra¬ 
tion. For instance, derivatives have been defined through 
a Fourier transformation to momentum space [a, m 1 
which is equivalent to the assumption that the functions 
on the mesh can be developed into a set of plane waves. 
In this method, the derivatives are quasi-exact for a given 
resolution of the mesh, and first- and second-order deriva¬ 
tives internally consistent. Similar ideas have been de¬ 
veloped in quantum chemistry under the label of discrete 
variable representation (DVR) [25l - l^ . A similar for¬ 
malism that provides an internal consistent scheme for 
derivatives and integrals is the Lagrange-mesh method 
that we will sketch in the following section. 


B. Lagrange-mesh representation 


The idea underlying the Lagrange-mesh method is that 
for each Gauss quadrature one can construct a set of 
basis functions for which orthogonality and completeness 
relations are exactly fulfilled when evaluated with the 
given quadrature [ifl, [H, [3 ■ This additional condition 
makes the LM method a srccial case of the slightly less 
rigorous concept of DVR [3 HI] ■ 

Lagrange meshes have been constructed for a multi¬ 
tude of different geometries and used for a wide range 
of applications, see (3 and references therein. We will 
use here the case of an equidistant 3d cartesian mesh. Its 
three directions are separable in the formalism, such that 
the presentation of the principles of the method for one 
dimension is sufficient. 

The underlying basis of a one-dimensional Cartesian 
equidistant Lagrange mesh is constructed as the set of 
functions Lpk{x) whose orthogonality relations are ex¬ 
act when evaluated with a simple 2A-point rectangular 
quadrature rule, sometimes called midpoint rule, [3 


dx (filix) (pk'(x) 

> dx'^ifilixr) ifk'ixr) = Skk' 


( 12 ) 


where dx is the distance between the collocation points 
located at 


Xr = r dx = ±dx/2, ±3 dx/2,..., ±(A — 1) dx/2 , (13) 

and where a and b are the boundaries of the numerical 
box [a, b] = [—Ndx, +Ndx]. A convenient representation 
of the 2N basis functions ipk{x) are plane waves of the 
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form 

<^>k{x) = x) , (14) 

where L = 2Ndx is the length of the numerical box and 
where A: = ±i, ±|, The real part of the 

ipk{x) is symmetric, has nodes on the boundaries of the 
box and a maximum at the origin, whereas their imag¬ 
inary part is skew-symmetric, consequently has a node 
at the origin, and maxima on the boundaries of the box. 
This also implies that (pl(x) = (p-k(x). 

The (fk(xr) form a complete set of functions to describe 
any function on the mesh points 

dx'^ifKXr) = Srs • (15) 

k 

Note that the box size L is not a multiple of the wave¬ 
length of the basis functions. Instead, twice the box size 
is an odd multiple of the wavelengths that take the val¬ 
ues 2L = 2L/1, 2L/3, 2L/5, ... 2L/{2N - 1). Both the 
real and imaginary parts of all plane waves in Eq. m 
are non-zero at all mesh points. 

As recalled in Ref. [^, in cartesian DVR and LM 
coordinate-space methods where the derivatives are de¬ 
fined through an expansion in plane waves, the analy¬ 
sis of a calculation’s infrared and ultraviolet cutoffs in¬ 
troduced by the basis is straightforward. This has to 
be contrasted with the much more involved analyses re¬ 
quired when working with an HO basis It has 

also been argued in Ref. [1^ that a DVR or LM repre¬ 
sentation of the nuclear many-body problem covers the 
relevant part of the phase space with a much smaller 
number of basis states than required by an HO basis. 


In practice, however, HO bases typically used for self- 
consistent mean-field calculations are much smaller than 
the typical number of mesh points used in the same kind 
of calculation. For a box with 20 points in every direction 
the number of linearly independent states is 64000, to be 
compared with a harmonic oscillator expansion with 20 
shells, which contains I4I68 states. 

While the basis functions ipk{x) of Eq. (fH!) are use¬ 
ful to discuss the mathematical properties of the LM 
method, the actual coordinate representation then em¬ 
ploys the set of 2N La gra nge interpolation functions 
fi{x) obtained as (l^. 1^.1^ 

fr{x) = dx'^(pl{Xr) ipk{x) 
k 

^ 1 sin[^(a:-rdx)] 

By construction, the Lagrange interpolation functions 
have the property to be equal to one at the mesh point 
Xr = rdx, and zero at all others, frixg) = i5r..|10l.[^.[^. 
When developed into the Lagrange functions, any func¬ 
tion on the mesh 

='^(l){Xr) fr{x) ='^(l>r fr{x) (17) 

r r 

is then simply represented by its values 4>r = ‘Pixr) at 
the 2N mesh points. 

The Lagrange functions are smooth and infinitely 
derivable. They can be used to define matrices repre¬ 
senting the first and second derivatives of functions dis¬ 
cretized througlS Eq. (fT7)l 


„(i) ^ df{x) 
dx 


{-ly- 
0 

^( 2 ) ^ d'^Mx) 

“ dx^ 


(2N)dx sin(7r(i — j)/(2A^)) 


for i^j, 
for i = j, 


(_I).-.+i 2 ( ^ ' cos[7r(z-j)/(2jV)] ^ . 

V(2IV)da:7 sin^ [Tr{i - j)/{2N)] 

for i = j. 


Mx^ 


1 - 


{2Ny 


(18) 


(19) 


The first derivative of any function (l){x) on the mesh is 

^ Unfortunately, the corrections of these expressions as given in obtained by multiplying the 2N X 2N matrix Drs by the 
the corrigendum to Ref. fT3l still contain a typographical error: 
the formula for the second derivative has a superfluous factor of 
two when i ^ j. 
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vector (jjr 

= = ( 20 ) 
S 

and similar for the second derivatives. Note that the 
derivative matrices have the property 
by construction [2^, which is not the case for hnite- 
difference formulas. As the derivatives of Eqs. (HI and 
(HU) correspond to full 2N x 2N matrices, their applica¬ 
tion is more time consuming than hnite-difference deriva¬ 
tives that correspond to a sparse band matrix. 

The full cartesian 3d representation of a function $(r) 
is then provided by 

( 21 ) 

pqs 

where the number of discretization points does not have 
to be the same in each direction. In that case, the deriva¬ 
tive matrices in Eqs. m and (HU) have to be set-up 
separately for each direction. 

As pointed out in Refs. a variational calcula¬ 

tion using a DVR or LM derivatives delivers very precise 
values for the total energy in spite of the individual ma¬ 
trix elements being much less accurate. In what follows, 
we will illustrate that this property implies very accurate 
total energies while separate terms of the Skyrme EDE 
are less well represented. In addition, we will show that 
using a LM results in a variational calculation. 

IV. NUMERICAL CONSIDERATIONS 

A. Numerical parameters of parameterizations 

Unless explicitly stated, we have used the SLy4 param¬ 
eterization. To explore the dependence of the numerical 
accuracy on the EDE, we have in addition tested a rep¬ 
resentative set of Skyrme parameterizations, as listed in 
Appendix [A] 

In the next two sections, we present calculations for 
doubly-inagic spherical nuclei ^°Ca, ^^^Sn and ^°®Pb, the 
neutron-rich nucleus ^^Ne, Cd, Sn and Te isotope chains 
and the fission path of ^"^°Pu. It is worth noting that 
we only included pairing for the isotopic chains and for 
^"^^Pu, see appendix m for details. In all other cases, 
pairing has been neglected. 

In appendix [C] we comment on the precise physical 
constants used during our calculations. 

B. Measuring accuracy 

The accuracy of a coordinate space calculation is lim¬ 
ited by the size of the box, the discretization length dx 
and the way derivatives and integrals are calculated. In 
order to properly judge these effects we employ two ways 
of analyzing results. Eor spherical nuclei we can compare 



FIG. 1. (Color online) Comparison between the errors on the 
total energy of obtained in calculations using different 

combinations of formulas for the derivatives (see text). The 
differences are taken with respect to the total energy obtained 
with Lenteur [2ll] . 

our 3d results with a one-dimensional spherical code that 
also represents the single-particle wave functions in coor¬ 
dinate space. Because of spherical symmetry, we can use 
extremely fine discretizations and the results can thus be 
considered exact to very high precision. Eor this purpose 
we use Lenteur [2l[ as a reference. 

Eor deformed nuclei, we no longer have access to such 
a comparison. Here we have to resort to looking at 3d re¬ 
sults as a function of both box size and mesh spacing: we 
compare results in small boxes with large mesh spacing 
to results in very large boxes with very fine mesh spacing. 


C. The use of derivatives and the variational 
principle 

The numerical cost of using LM derivatives is much 
higher than the ED alternative. To control the compu¬ 
tational time, three options have been considered: they 
differ by the way derivatives are calculated during the 
mean-held iterations and after convergence. The hrst 
option (ED-I-ED) has been used in the hrst applications 
of the codes Q where derivatives were exclusively cal¬ 
culated by ED. The second one (FD-f-LM) is the most 
used one for more than 20 years: ED derivatives are used 
during the iterations but the energies are recalculated 
after convergence by LM formulas. Finally, in the last 
option (LM-I-LM), the LM formulas are used during the 
iterations and after convergence. 

In practice, we use a seven-point difference formula for 
the hrst order and a nine-point formula for the second or¬ 
der derivatives when employing ED formulas. It has been 
shown earlier in Ref. [34| that this provides an efficient 
compromise in terms of overall speed and precision. 

Figure[T]illustrates the accuracy on the total energy ob- 
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Energy (MeV) 

dx = 1.0 fm 

dx — 0.83 fm 

dx = 0.71 fm 

dx = 0.549 fm 

Lenteur 

Kinetic -f c.m. correction 

3908.548 

3880.647 

3872.035 

3867.506 

3866.190 

Direct Coulomb 

831.801 

829.241 

828.433 

828.004 

827.876 

Coulomb Exchange 

-31.415 

-31.319 

-31.289 

-31.273 

-31.269 

Ep2 

-22749.578 

-22510.048 

-22435.391 

-22395.941 

-22384.379 

Epr 

1368.924 

1343.506 

1335.611 

1331.436 

1330.206 

Ep2 + c 

14812.851 

14631.848 

14575.413 

14545.584 

14536.832 

Ep^p 

323.771 

318.147 

316.431 

315.539 

315.287 

Epvj 

-99.730 

-97.595 

-96.917 

-96.554 

-96.445 

Total Skyrme Energy 

-6343.762 

-6314.142 

-6304.854 

-6299.937 

-6298.501 

Total Energy 

-1634.828 

-1635.574 

-1635.675 

-1635.700 

-1635.703 


TABLE 1. Decomposition of the total energy between the terms of the Skyrme parametrization SLy4 for using the 

FD+LM option. All energies are in MeV. See Sect. Ill Al for the definition of the various terms. Values obtained with the 
spherical Id code Lenteur are given for comparison. 


Energy (MeV) 

dx = 1.0 fm 

dx = 0.83 fm 

dx = 0.71 fm 

dx = 0.549 fm 

Lenteur 

Kinetic -f c.m. correction 

3866.323 

3866.165 

3866.182 

3866.182 

3866.190 

Direct Coulomb 

827.922 

827.889 

827.882 

827.878 

827.876 

Coulomb Exchange 

-31.269 

-31.268 

-31.268 

-31.268 

-31.269 

Ep2 

-22384.936 

-22384.188 

-22384.322 

-22384.320 

-22384.379 

Epr 

1330.300 

1330.193 

1330.201 

1330.200 

1330.206 

Ep2 + a 

14537.174 

14536.691 

14536.789 

14536.787 

14536.832 

Ep^p 

315.238 

315.275 

315.284 

315.284 

315.287 

Epvj 

-96.430 

-96.444 

-96.445 

-96.445 

-96.445 

Total Skyrme Energy 

-6298.657 

-6298.473 

-6298.493 

-6298.494 

-6298.501 

Total Energy 

-1635.678 

-1635.687 

-1635.696 

-1635.700 

-1635.703 


TABLE 11. Same decomposition as in Table [II but for the LM+LM option. 


tained using these three options. The LM + LM choice 
is by far the most accurate. As it can be seen in Ta¬ 
ble m the result obtained with a mesh size of 1.0 fm 
differs by only 25 keV from the Lenteur result. The 
FD-fLM option is less accurate, but already sufficient 
for most applications, with an error of around 100 keV 
for dx = 0.8 fm. It is better by nearly an order of magni¬ 
tude than the FD-t-FD choice. Results presented in the 
following have been obtained with the FD-f-LM option, 
except otherwise stated. 

Both the FD-I-LM and LM-I-LM calculations underesti¬ 
mate the binding energy, as it should be for a variational 
calculation. This is due to the fact that the single-particle 
wave functions are expanded on a complete and closed 
basis for given box size and mesh discretization length, 
see Eq. (HI. Increasing the box size and/or decreasing 
the mesh discretization length enlarge the accessible sub¬ 
space of the Hilbert space and lead to a monotonous 
convergence of the energy. By contrast, such a basis can¬ 
not be defined for the FD-I-FD option for which the cal¬ 
culation systematically overestimates the binding energy 
of 208pb. 


The same applies to mesh calculations with Fourier 
derivatives, as can be deduced from the convergence anal¬ 
yses in Refs. sm. While for a given dx the overall ac¬ 
curacy of the binding energy found there is very similar 
to the one we find for LM-I-LM calculations, the energy 
does not converge monotonically when decreasing dx. 

While the use of LM derivatives after having used FD 
ones during the iterations (FD-fLM) is sufficient to ob¬ 
tain an upper bound of the total energy since any wave 
function discretized on a mesh can be expanded on the 
LM basis, the errors on the various individual terms of 
the Skyrme EDF can be very large, as can be seen in 
Table U While the total energy varies by slightly less 
than one MeV when dx is decreased from 1.0 fm to 
0.549 fm, the variation in the kinetic energy is of the order 
of 40 MeV, counterbalanced by a similar change in the 
Skyrme energy. The situation for the LM-fLM scheme 
is shown in Table m It indicates a similar effect, but on 
a much smaller scale: the total energy varies by 20 keV 
while the kinetic energy varies by roughly 150 keV. 

When performing symmetry restoration and configu¬ 
ration mixing by the GCM, a high level of accuracy is 
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FIG. 2. (Color online) Energy difference between a reference 
calculation performed with 23 points and calculations per¬ 
formed with N points for the three nuclei ^°Ca, ^®^Sn and 
^°®Pb. In all cases, the step size is equal to 1.0 fm. 


required to avoid buildup numerical noise while solving 
the Hill-Wheeler-Griffin equation. This calls for the use 
of LM derivatives in these calculations, as done since our 
first applications [3^ . 


D. Determining box sizes and mesh spacings 

The first requirement of a coordinate space calculation 
is that the box in which the nucleus is confined is large 
enough to avoid any spurious effect due to the truncation 
of the wave functions. The influence of the box size on 
the total energy for three spherical nuclei is represented 
in Fig. [21 The same mesh size da; = 1.0 fm is used in all 
calculations while the number of discretization points is 
varied, changing thus the volume of the box. The cal¬ 
culation in the largest box, using 23 points, is taken as 
a reference. The errors decrease quickly when the box 
size is enlarged. If one requires that the error is smaller 
than a keV, we see that taking boxes with half-sides of 
11 fm for ^°Ca, 15 fm for ^^^Sn and 20 fm for ^°®Pb is 
sufficient. Since the numerical effort required for ^°Ca is 
very low we opted to use a slightly larger half-side of 13 
fm in order to further increase our accuracy to about 0.1 
keV. 

Similar analyses have been performed for all nuclei con¬ 
sidered in this paper. Since several nuclei in the isotopic 
chains around Z = 50 are deformed, we have performed 
all calculations with the same box size as This 

choice allows us to calculate all isotopes with the same 
numerical conditions. The box dimensions are summa¬ 
rized in Table M The columns Cx , Cy and Cz indicate 
the size of the box in which the Coulomb problem is 
solved. For every system, the box size was varied for 
fixed dx until the energy did not change by more than 
0.1 keV, with the exception of the ^^°Pu for which this 
limit was of 1 keV. 


nuclei 

II 

Lz 

II 

Cz 

«Ga 

26 

26 

26 

26 


30.8 

30.8 

46.8 

46.8 

208pb 

40 

40 

60 

60 

Z « 50 

40 

40 

60 

60 

240p^ 

40 

60 

80 

120 


TABLE III. Edge lengths in fm of the boxes used to solve the 
self-consistent mean-field equations {Ly) and to determine the 
Coulomb potential {Cy) for the nuclei studied in this paper. 
Depending on the symmetries imposed on the nucleus, only 
half of of the length is treated numerically in most cases. 



FIG. 3. (Color online) Error on the total energy using the 
FD+LM option as a function of the number of mean-field 
iterations for “^^Ca for different values of dx. Calculations were 
initialized with Nilsson-model single-particle wave functions. 
The box length is 26 fm. 

A non-ambiguous comparison between calculations 
performed with different mesh discretizations dx can only 
be achieved when the volume of the box is conserved. 
This is realized by determining the value of dx in such a 
way that the box has the same size for each number of 
mesh points. 


E. Convergence of the iterative procedure 

Decreasing the mesh size improves the accuracy. How¬ 
ever, this has a price in computing time. First, to keep 
the same box size requires to increase the number of dis¬ 
cretization points. A second factor increasing the com¬ 
puting time is that the time step of the imaginary-time- 
step method [H, implemented in the codes [T^ has 
to be decreased with decreasing mesh size. This consid¬ 
erably slows down the convergence. In Fig. O we show 
the evolution of the error on the total energy relative to 
Lenteur during the iterations for the nucleus ^'^Ca for 
different mesh discretizations dx. The most accurate re¬ 
sult after 100 iterations is obtained with dx = 1.0 fm. 























FIG. 4. (Color online) Differences between the total energy of 
calculated on a 3d mesh with different approximations 
in the calculation of the Coulomb energy and that obtained 
with Lenteur. The three lines correspond to using 3-, 5- and 
7-point FD formulas for the calculation of the Laplacian in 
the Poisson equation. 

Gaining an order of magnitude of accuracy after conver¬ 
gence requires to carry out roughly 100 more iterations 
for the step sizes represented on the figure. 

F. Treatment of the long-range Coulomb 
interaction 

The direct Coulomb energy requires a special treat¬ 
ment because of its long range. One of the spatial in¬ 
tegrations in Eq. 0 can be eliminated through the cal¬ 
culation of the Coulomb potential of the protons, which 
satisfies the electrostatic Poisson equation 

AU{r) = -47re^pp(r), (22) 

where is the square of the elementary charge. When 
solving this equation, boundary conditions need to be 
imposed at the edge of the box.These can be easily 
constructed when recalling that at large distances the 
potential is entirely determined by the multipoles of 
the nucleus’ charge distribution {Qim)- Expanding the 
Coulomb potential on spherical harmonics and keeping 
terms up to £ = 2, the Coulomb potential outside the 
box is approximated by 

I ^2 (Q2o)h2o(r) + (O22) ^6^22(1') ^ 23 ) 

which provides the boundary condition for the numerical 
solution of Eq. (I22p . The direct Coulomb energy is then 
calculated as 

Econ\ = \ Jd^rU (r) pp{r). (24) 


As for the nuclear part of the energy, the accuracy of the 
electrostatic potential, obtained by solving Eq. (l22ll . is 
limited by three factors: the size of the box, the mesh 
discretization length dx and the way derivatives are cal¬ 
culated. 

A suitable box size for the Coulomb problem has to be 
larger than for the Skyrme EDF. This is a direct conse¬ 
quence of the long range of the Coulomb force. To make 
negligible the contributions to the boundary conditions 
of terms higher than £ = 2, see Eq. (l23ll . one has to calcu¬ 
late the Coulomb potential in a box larger than the one 
used for the nuclear part of the interaction. Typical val¬ 
ues are given in Table Hill For light nuclei such as ^°Ca, 
no extra points for Coulomb need to be added, while the 
box has to be significantly enlarged for heavier systems 
in the ^^^Sn and ^*^®Pb regions. For the calculation of the 
fission barrier of heavy nuclei such as ^^*^Pu up to very 
large deformations, the Coulomb box size has to be two 
times larger than the one needed for the Skyrme EDF to 
obtain the same nuclear accuracy on all the energies. 

The Laplacian in Eq. (1221) has to be approximated on 
the mesh in such a way that the accuracy on the Coulomb 
energy is similar to the one of the other terms in the EDF. 
We show in Fig. 2] the gain in accuracy on the total en¬ 
ergy of ^°®Pb obtained by going from a three-point to a 
seven-point FD formula for the Laplacian. Already a hve- 
point formula brings the required accuracy and is used in 
all other calculations reported here. One can easily un¬ 
derstand that a lower-order finite-difference formula than 
the one used to calculate the kinetic energy is sufficient 
for the Laplacian in Eq. (l22l) : the typical length scale 
of the variation of the Coulomb potential is much larger 
than the scale on which the wave functions vary. 

The final factor for the accuracy of the Coulomb solu¬ 
tion is the mesh discretization length dx. As the effect 
of the Coulomb term is already incorporated in all of the 
applications, we will not discuss it separately. 


V. DISCUSSION 
A. Binding energies 

Provided that the box size is large enough, the main 
factor determining the accuracy of our implementation of 
mesh calculations is the discretization length. In Fig. O 
the energy difference with respect to Lenteur results is 
plotted for three doubly-magic spherical nuclei, as a func¬ 
tion of the mesh discretization dx for a representative set 
of Skyrme parameterizations. It is remarkable that the 
interactions are grouped according to their effective mass 
(see Appendix [A] for the actual values): interactions with 
larger effective mass m* give systematically more accu¬ 
rate results than interactions with smaller ones. This 
property is related to the term Epj- term of the Skyrme 
EDF in Eq. (I2|) that, in our experience, is the least well 
represented on a mesh. Since the magnitude of this term 
increases when the effective mass decreases, the accuracy 
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dx (fm) 


FIG. 5. (Color online) Differences between the total energy calculated with our 3-dimensional code and with the spherical code 
Lenteur for “^^Ca (a), ^®^Sn (b) and ^°®Pb (c), as a function of mesh distance dx. Results are plotted for a representative set 
of Skyrme parameterizations, without pairing. Results for SLy5, T22 and T65 are not shown, but are indistinguishable from 
the SLy4 results on the scale of this graph. 


obtained for a given mesh size deteriorates for lower ef¬ 
fective mass. 

One can see that the accuracy obtained with a mesh 
discretization as large as dx = 1.0 fm is lower than 
1.0 MeV for ^°®Pb. The energy difference decreases 
to a few hundred keV for dx = 0.8 fm and to a few 
keV for dx = 0.6 fm. Note that a similar accuracy for 
dx = 0.6 fm was found for a 2-dimensional code based on 
splines To obtain an agreement between the spheri¬ 
cal code Lenteur and our 3-dimensional codes below the 
1 keV level would require to increase the box size but 
also to make the codes more similar. For a nucleus with 
a binding energy larger than 1 GeV, this implies a rela¬ 
tive discrepancy better than 10“^ and there are several 
sources of differences in the codes that can play a role, 
none of which is easy to control. 


B. Deformation energy curves 

Let us now study the convergence properties of our 
numerical scheme for the fission path of Our mo¬ 

tivation is twofold: ^"‘^Pu is a frequent benchmark for 
models that describe fission [sol - El l but also for numeri¬ 
cal algorithms Sin. The energy curve of this nucleus 
presents two minima at prolate deformations, the ground 
state and a fission isomer. In Fig. [6l we show the varia¬ 
tion of the energy with deformation. 

The box used for these calculations has the same size 
for all discretizations, as indicated by Table IIIII When 
the left-right symmetry is broken, the number of points 
along the z direction is doubled. We have performed 
calculations with four different mesh discretization, dx = 
1.0, 0.82, 0.69 and 0.60 fm and tested the convergence as 
a function of dx by taking the difference with respect to 
the results obtained with dx = 0.6 fm. For each value of 
dx, the energy at each deformation is the energy relative 
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FIG. 6. (Color online) Energy curve for ^^’’Pu calculated with 
dx = 0.6 fm. The regions where the deformation is axial, 
triaxial or axial and reflection asymmetric are indicated in 
blue, green and red respectively. 


to the prolate ground state. 

The energy curve obtained with dx = 0.6 fm is shown 
in Fig.[6l The topography obtained for other values of dx 
is the same. Shapes are triaxial in the vicinity of the first 
barrier, whereas everywhere else they remain axial. At 
deformations smaller than the one of the fission isomer 
the configurations are reflection symmetric, whereas at 
larger deformations they are increasingly asymmetric. 

We will use this curve as a reference to determine the 
accuracy of the calculations carried out for other values 
of dx. For each dx, the ground state energy is taken as 
the zero of the energy. The results are shown Fig.[71 The 
properties of the minimum are summarized in Table El 
The error decreases roughly by an order of magnitude by 
going from dx = 1.0 to 0.82 fm and from dx = 0.82 to 




















10 



FIG. 7. (Color online) Energy differences between the results 
obtained for with dx = 1.0, 0.82 and 0.69 fm and those 
corresponding to dx = 0.6 fm. 



FIG. 8. (Color online) Absolute difference in neutron rms 
radius for different Skyrme parameterizations for ^^Ne for 
dx = 0.8 fm. The reference calculation was performed in 
a box with N = 20. 


0.69 fm. At dx = 1.0 fm the error is of the order of a few 
times 100 keV, with a rather large oscillation. For a mesh 
discretization of 0.82 fm, the error becomes lower than 
100 keV (except in the vicinity of the spherical configu¬ 
ration where it reaches 150 keV, but this configuration is 
very excited) and is quite acceptable for the calculation 
of energy curves. Decreasing still the discretization to 
0.69 fm reduces the error to values around a few times 
10 keV at most. 

Some published results allow for a comparison between 
the accuracies of mesh calculations and of calculations 
using an expansion on an HO basis. Pei et al. Q have 
performed calculations on an axial mesh using B-splines 
and on HO bases either spherical or deformed, with 20 
oscillator shells in both cases. The accuracy obtained 
in 0 on a mesh of dx = 0.65 fm seems very similar to 
the one we obtain. The use of a spherical HO basis is 
rather unreliable, with an error larger than 1 MeV al¬ 
ready for the excitation energy of the fission isomer and 
that quickly increases to several MeV at larger deforma¬ 
tions. For an axial oscillator basis, the results are similar 
to those that we obtain with a mesh size of 0.82 fm up 
to the first barrier but the accuracy deteriorates rapidly 
for larger deformations, being of several hundreds of keV 
at the deformation corresponding to the fission isomer. 
Similar results can be found in for ^®"^Hg and in 
for 256Fm. 

As a number of shells significantly larger than 20 is 
numerically prohibitive, one either has to resort to a 
two-center oscillator basis or one has to construct a suit¬ 
able subspace within a much larger one-center HO basis 
by carefully selecting the low-lying single-particle states. 
The former option is developed in Ref. whereas the 
latter has been used during the construction of the unedfl 
parametrization [i^, where the lowest 1771 basis states 
out of a basis of 50 HO shells has been kept. The accu¬ 
racy obtained in this way for the excitation energy of the 


E 

VI— 



Radius (fm) 


FIG. 9. (Golor online) Radial density profile of ®^Ne in dif¬ 
ferent box sizes with dx = 0.8 fm. 


fission isomer is of the order of 100 keV. As a comparison, 
the experimental excitation energy of the fission isomer 
that can be found in the literature is 2.4 ± 0.3 MeV M, 
approximately 2.8 MeV [5l| and 2.25 ± 0.20 MeV sj. 
In the light of these error bars, a numerical accuracy of 
100 keV is sufficient for the adjustment of an EDF. How¬ 
ever, from the published results by Pei et al. Q, it can be 
estimated that the numerical error on the fission barrier 
height is a few times these 100 keV. Similar results have 
been obtained in the case of the RMF method [l^, • 


C. Radial density distribution 

The rms radius is intimately linked to the radial den¬ 
sity distribution of a nucleus. One can expect that it 
is particularly sensitive to the box size for nuclei with 
a large excess of neutrons. Tests have been performed 
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dx (fm) 


FIG. 10. (Color online) Absolute difference between the total rms radii calculated on a 3d mesh with respect to those of 
Lenteur as a function of the step size da:for the spherical nuclei (a), ^®^Sn (b), and ^'^®Pb (c) and calculated with the 

Skyrme parameterizations as indicated. 
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FIG. 11. (Color online) Isotopic shifts 5r^{N,Z) with respect 
to ^^^^Sn for different Sn isotopes and different mesh sizes. 


for the very neutron-rich nucleus ^'^Ne by varying the 
box size for a fixed mesh discretization dx = 0.8 fm. To 
avoid any ambiguity in the calculation, pairing has been 
omitted. The results are presented in Fig. [8l where we 
show the difference in total rms radius as a function of 
the box size for a representative set of EDF parameter¬ 
izations. For the size of the box recommended for ^°Ca 
in Table [TTTl the number of points is 16 for a mesh size 
of 0.8 fm. It leads to an error of the order of 10“^ fm for 
most interactions, the results being slightly less accurate 
for SV-min. For smaller boxes, the accuracy of radii is 
lower and depends on the interaction. 

In Fig.|9]the radial profile of the total density of ^^Ne is 
plotted as a function of the box size. The distortion of the 
density in the smallest box is large and demonstrates that 
half the box size has to be larger than 8.0 fm. In all other 
boxes, the exponential tail of the density distribution is 
well described, up to the point before the last one. For 
a box size around 12 fm, the density is well described 


up to a decrease of the central density by six order of 
magnitudes. 

The confinement in a volume is less evident in an ex¬ 
pansion on a basis than in a mesh calculation, but it is 
also present. While oscillator basis functions extend to 
infinity, they are in practice strongly localized by their 
Gaussian form factor. If one takes its classical turning 
point as a measure for the extension of a HO state, one 
obtains, for ^°®Pb and 20 oscillator shells, a value for the 
turning point that varies from 14 fm for £ = 0 to 16 fm 
for i = 20h. To increase the value of this turning point to 
20 fm would require to use 28 oscillator shells for £ = 0. 
This effect of confinement by an oscillator basis has been 
put in evidence in Ref. for the case for ^^^Zr. 

For comparison, the experimental uncertainty on rms 
charge radii for the Ne isotopes (up to A = 28) varies 
from 0.002 fm close to stability to 0.02 fm for exotic iso¬ 
topes [HI. It is interesting to note that the numerical 
accuracy of a mesh mean-field calculation has a similar 
level (provided the box is large enough), but that the 
model already introduces uncertainties on the rms radii 
that are at least one order of magnitude larger Q. 

In Fig. [TOl we compare the total rms radii calcu¬ 
lated with decreasing mesh sizes to those obtained with 
Lenteur for three spherical nuclei "^^Ca, ^^^Sn and ^°®Pb. 
The agreement is already very satisfying for a large mesh 
size of 1.0 fm, with one order of magnitude gained in ac¬ 
curacy when decreasing the mesh size to 0.8 fm, which is 
the usual value of production calculations. An interest¬ 
ing feature that cannot be deduced from Fig. [T0]is that 
all of the parameterizations, with the exception of un- 
edfO, always produce an rms radius that is smaller than 
the Lenteur result. 

In Fig. [TTl we present the isotopic shifts Sr'^{N, Z) for 
a range of even-even Sn nuclei, the reference being ^^^Sn. 
All curves almost exactly coincide. This demonstrates 
that the isotopic shifts are quite reliable even with coarse 
meshes. Similar results are obtained for Cd, Xe and Te 
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isotopes. 


D. Two-neutron separation energies 

To put into evidence changes of nuclear structure with 
nucleon number, one often uses mass filters that are com¬ 
puted by taking specific differences between the binding 
energies of neighboring nuclei. The simplest filter is the 
two-nucleon separation energy that is defined as the en¬ 
ergy difference between two isotopes (or isotones) whose 
nucleon number differs by two. In Fig. 1121 we show the 
evolution of the two-neutron separation energies, S 2 m 
of even-even nuclei for three neighboring isotopic chains 
when the mesh size dx is decreased. For each discretiza¬ 
tion dx we have plotted the difference of the S 2 n values 
to the one obtained at dx = 0.63 fm. Even with a mesh 
size as large as dx = 1.0 fm, the accuracy on the S 2 n is 
already better than 100 keV, which is small enough for 
most applications. The mesh size used in most of our 
published applications, dx = 0.8 fm leads to an accu¬ 
racy better than 10 keV. In the bottom panel, the two- 
neutron separation energies of the three isotope chains 
are plotted for four values of dx. The curves cannot be 
distinguished using a scale adapted to the variation of 
S 2 n as a function of the neutron number. This result 
is in strong contrast with respect to some published cal¬ 
culations using an expansion on an oscillator basis [s^, 
where special algorithms have to be devised to smooth 
numerical irregularities that can be of the order of few 
hundred keV. 


E. Multipole moments 

The dimensionless ground-state quadrupole moments 
/32 of even-even Te isotopes are shown in Fig. [131 Dif¬ 
ferences between the curves corresponding to different 
values of dx are tiny and not signihcant. Similar results 
were obtained for the Cd and Sn isotopes. 

We now examine how the multipole moments of ^^°Pu 
along the fission path are affected by the mesh size. In 
Figs. [HI and [T51 we show the octupole and hexadecapole 
moments, respectively, in the region of the fission path 
where parity is broken. Similar results obtained for the 
axial and triaxial cases are not shown. In Tables IIVI 
and El we show the multipole moments of the ground 
state and fission isomer of ^^°Pu for the different mesh 
discretizations as obtained by unconstrained calculations. 

From Figs. [HI and [T5] we see that the overall sequence 
of shapes along the fission path is robust with respect to 
the mesh spacing. The fission path is already precisely 
defined at the coarsest mesh {dx = 1.0 fm) we used. 
A single exception can be seen at the onset of octupole 
deformation; in the vicinity of this point, however, the 
energy surface is very flat in direction. 

On a smaller scale, the multipole moments do vary as 
a function of the mesh discretization. This is best visible 
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FIG. 12. (Color online) Absolute differences of two-neutron 
separation energies between four mesh discretizations and a 
calculation with a mesh size of dx = 0.63 fm for Cd (a), Sn (b), 
and Te (c) isotopes. In (d) we plot the two-neutron separation 
energies of all of these isotope chains as a reference. 
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dx (fm) 

E (MeV) 

P 2 

Pi 

Pe 

Ps 

Pio 

1.0 

-1801.909 

0.288 

0.160 

0.043 

-0.002 

-0.003 

0.849 

-1802.770 

0.289 

0.163 

0.045 

0.001 

0.001 

0.739 

-1802.929 

0.292 

0.164 0.046 

0.002 

0.002 

0.653 

-1802.969 

0.290 

0.165 

0.046 

0.002 

0.001 


TABLE IV. Properties of the ground state of as ob¬ 

tained from non-constrained calculations. 


dx (fm) E (MeV) /32 134 Pe Ps Pio 


1.0 -1796.929 

0.849 -1797.950 
0.739 -1798.099 
0.653 -1798.123 


0.832 0.494 0.344 0.279 0.255 
0.840 0.510 0.367 0.303 0.278 
0.847 0.528 0.388 0.319 0.268 
0.841 0.516 0.375 0.312 0.259 


FIG. 13. (Color online) Mass P 2 quadrupole moment as a 

function of neutron number for a series of Te isotopes. TABLE V. Properties of the superdeformed fission isomer of 

as obtained from non-constrained calculations. 
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FIG. 14. (Color online) Mass octupole moment P 3 along the 
fission path for ^'‘“Pu. 



FIG. 15. (Color online) Mass hexadecapole moment P 4 for 
the parity breaking configurations along part of the fission 
path for ^'‘'’Pu. 


FIG. 16. (Color online) Eigenvalues of the single-particle 
Hamiltonian for the ground state of ^“^^Pu as a function of 
mesh discretization dx. Only neutron single-particle levels 
within 1.5 MeV of the Fermi energy are shown. 

in Tables HVl and fVl Since our method hinges on the vari¬ 
ation of the total energy in Eq. o, there is no guarantee 
that the values of the multipole moments converge in a 
predictable way. It is, however, reassuring to see that 
the typical variation of these moments is of the order of 
a few percent to at most about 10 percent. The larger 
variations present themselves in the higher-order Pg,Ps 
and Pio moments. These are more difficult to resolve on 
coarse meshes because of the high number of nodes their 
associated Legendre polynomials have. 


F. Single-particle levels 

In Fig. [in] we show the evolution of the neutron single¬ 
particle levels within 1.5 MeV of the Fermi energy in the 
ground state of ^^°Pu as a function of the mesh spacing 
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dx. While slight shifts of the position of the levels are 
observed as a function of the mesh size, the largest error 
at da; = 1.0 fm is of the order of 100 keV. One can also 
note that the level ordering within the parity subspaces 
is the same for all values of dx. A similar dependence 
on box parameters is found for the proton states and the 
lighter nuclei studied here. 


VI. CONCLUSION 

The aim of this paper was to study the numerical accu¬ 
racy of the solution of the self-consistent mean-field equa¬ 
tions using a discretization on a 3-dimensional cartesian 
coordinate-space mesh. Three elements permit to control 
its numerical accuracy. The first one is the method used 
to calculate derivatives. Using Lagrange-mesh deriva¬ 
tives leads to much more accurate results than hnite- 
difference formulas for derivatives. In addition, a carte¬ 
sian Lagrange mesh corresponds to a representation in a 
closed subspace of the Hilbert space, such that it always 
provides an upper bound to the binding energy that be¬ 
comes tighter when adding points outside a given box or 
when decreasing the distance of mesh points in a given 
box. Neither is the case for hnite-difference derivatives. 
However, we have shown quantitatively that the accuracy 
of a calculation that uses finite-difference formulas during 
the iterations can be significantly improved upon by re¬ 
calculating the EDF at convergence with Lagrange-mesh 
derivatives. Again, this procedure provides us with an 
upper bound of the energy, thus restoring the variational 
character of the calculation. Using Lagrange derivatives 
during the iterations allows to still improve the accuracy 
on energies but at the cost of at least doubling the com¬ 
puting time. 

The second element on which mesh calculations de¬ 
pend is the size of the box in which the nucleus is con¬ 
fined. The examples of doubly-magic nuclei and neutron- 
rich ^"'Ne illustrate that results for energies and densities 
are already stable at small box sizes. 

Thirdly, the quality of the results depends on the mesh 
size with errors on energies that are almost independent 
on the number of neutrons and protons and on the shape 
of the nucleus. A mesh size dx = 0.8 fm guarantees an 
accuracy that is in general better than 100 keV, which 
corresponds to a relative accuracy of less than a tenth 
of percent, even for lighter nuclei. Decreasing the mesh 
size to 0.7 fm permits to gain nearly an order of magni¬ 
tude and to reach an accuracy that is well below all the 
uncertainties of the mean-field model. 

One can summarize these results by concluding that a 
mesh technique as implemented in our codes is flexible 
(it can accommodate any kind of symmetry breaking), 
robust (the accuracy can be controlled by an adequate 
choice of the three elements mentioned above) and that 
it can be very accurate if needed. The positive aspect of 
our numerical scheme is that using a mesh size of 0.8 fm, 
as used in most of our past applications ensures an accu¬ 


racy better than 100 keV on energies and reliable shape 
properties for nuclei of any mass. 

Our study has been focussed on the solution of the 
mean-field equations and we have not touched the de¬ 
scription of pairing correlations. There has already been 
a study of this problem by Terasaki et al. [^. It 
should be revisited today to take into account new de¬ 
velopments. However, the problem is not exclusively 
a problem related to the way the mean-field equations 
are solved. The description of single-particle states well 
above the Fermi energy is probably very different when 
using a discretization on a mesh or an expansion on an 
oscillator basis. 
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Appendix A: Parameterizations used 

We here give some remarks on the interactions used 
throughout the text. 

• SLy4 and SLy5 [s^: these parameterizations were 
used as intended. 

• T22 and T65 [s^: these parameterizations were 
used as intended. 

• UNEDFO [13 : while this parametrization was ad¬ 
justed with a non-zero pairing interaction, we used 
it without any pairing. 

• SV-min [6lj: this parametrization was adjusted 

with non-equal nucleon masses trip. Lenteur 

does not handle this option, so we used instead the 
average value as nucleon mass. 

• SkI3 [s^: this parameteriation was adjusted with 
the inclusion of a perturbative two-body c.m. cor¬ 
rection, an option not included in Lenteur. 

The values of their isoscalar effective mass are listed in 
Table O 
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Appendix B: Pairing interaction 


The density-dependent pairing interaction we used for 
the isotope chains and the fission path of is defined 
by M- 





1 — a 


po(R-) 

Ps 


d{r-r'), 


(Bl) 

where po{R) is the isoscalar density at R = i(r-l-r'). The 
parameters take the values a = 1, ps = 0.16 fm~^ and 
Vb = 1250MeVfm“^. In addition, this interaction was 
supplemented by two cutoffs, one above and the other 
below the Fermi energy, in order to eliminate the basis- 
size dependence of the total energy. They are defined by 
two Fermi functions 


(B2) 

where Ag is the Fermi energy, the single-particle energy 
of the single-particle state k and we chose pq = 0.5 MeV, 
and Aeg = 5.0 MeV for protons and neutrons. 


Appendix C: The role of physical constant when 
using Skyrme EDFs 


By default, the physical constants used in our calcula¬ 
tions are the following [s^: 

= 1.43996446 MeV fm, (Cl) 

m = = 938.9187125 MeVc"^ (C2) 

h^/{2m) = 20.735519104 MeV fm^ . (C3) 

Whenever possible, we used the value of that 

was used during the adjustment of the parametrization. 
It might seem superfluous to completely specify the phys¬ 
ical constant used, but the results of our calculations de¬ 
pend on the precise values of these constants. In particu¬ 
lar, the level of agreement between Ev8 and Lenteur de¬ 
scribed in Sect. IV Al is only attainable when these codes 
use exactly the same numerical values for the physical 
constants. 

In fact, significant errors can be introduced when the 
values of the physical constants are slightly changed. The 
seemingly innocuous value of j (2m) plays in fact a very 
important role. Figure [T3 shows Lenteur calculations for 
the spherical nuclei "‘^Ca, ^^^Sn and ^°®Pb with SLy4. 
Every point was calculated by slightly changing the value 
of fi^/(2m) from 20.73553 MeV fm^, the SLy4 value. We 
see that using a value for j{2m) that is not consistent 
with the value used during the fit of the EDF can lead to 
an error of several MeV on the total energy, k? /{2m) is 
after all the proportionality constant of the kinetic energy 
in Eq. O- Typical values that have been used over the 
years vary at least between 20.73 and 20.7363 MeV fm^. 
If the exact values of the physical constants used during 



FIG. 17. (Color online) Energy difference for the spherical 
nuclei ^®^Sn and for calculations with Lenteur 

using SLy4 with a modified value of f? /{2m). The reference 
calculation is that obtained using the value used during the 
adjustment of SLy4. 



FIG. 18. Difference between the total energy obtained with 
Lenteur for '**’Ca when using a rounded value for the density 
dependence parameter a in Eq. © instead of the full double 
precision value a = 1/6 of the SLy4 parameteriation. 


the adjustment of a given parametrization are not avail¬ 
able, then one cannot reliably compare the results with 
experimental data. In this case, one cannot judge the 
predictive power of this parametrization. 

Similar concerns arise for the parameters of the Skyrme 
interactions. The energy obtained in our calculations is 
more sensitive to some Skyrme parameters than to oth¬ 
ers, but the close agreement observed in Sect. IV Al is not 
obtainable without carefully checking that the Skyrme 
parameters are completely consistent across codes. That 
this is not trivial can be concluded from Fig.[THl There we 
plotted the relative difference in energy found by Lenteur 
between modified versions of the SLy4 functional and the 
correct SLy4. The interaction parameters are the same 
for every point, save for the density dependence param¬ 
eter a in Eq. There are only very few parameteriza- 
tions for which the value of a corresponds to a terminat¬ 
ing decimal, for example SV-min for which a = 0.255368. 
For the large majority of parameterizations the value of 
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Parametrization 

m*/m 

SLy4 & SLy5 

m 

0.68 

T22 & T65 


0.7 

SkI3 

[60] 

0.577 

SV-min 

m 

0.95 

unedfO 

m 

0.9 


TABLE VI. Value of the isoscalar effective mass in units of 
the nucleon mass m*/m of the interactions used throughout 
this paper. 


a is either 1/3 or, as in the case of SLy4, 1/6. Both of 


these correspond to a repeating decimal number, whose 
numerical representation might differ from code to code. 
Using a = 0.1667 in a calculation with SLy4 corresponds 
to a rounding error of a — 1/6 ~ 3.33 x 10“®, which in¬ 
troduces an error of the total binding energy of ^°Ca of 
a few tens of keV. 

It can clearly be seen that a limited representation of 
a implies a roundoff error that has a visible effect on 
the energy. This kind of error shows up when comparing 
Lenteur and Ev8 results and for this reason we conclude 
that relative errors smaller than 10“® become meaning¬ 
less. Similar analyses can be made for the other inter¬ 
action parameters, including the values of physical con¬ 
stants used to fit the interaction. 
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